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Abstract — This contribution considers one central aspect 
of experiment design in system identification. When a con- 
trol design is based on an estimated model, the achievable 
performance is related to the quality of the estimate. The 
degradation in control performance due to errors in the 
estimated model is measured by an application cost function. 
In order to use an optimization based input design method, 
a convex approximation of the set of models that satisfies 
the control specification is required. The standard approach 
is to use a quadratic approximation of the application cost 
function, where the main computational effort is to find the 
corresponding Hessian matrix. Our main contribution is an 
alternative approach for this problem, which uses the structure 
of the underlying optimal control problem to considerably 
reduce the computations needed to find the application set. 
This technique allows the use of applications oriented input 
design for MPC on much more complex plants. The approach 
is numerically evaluated on a distillation control problem. 

I. Introduction 

System identification for control concerns the problem of 
using experimental data from a dynamical system to identify 
a model to be used for control design, see e.g. [1], [2], 
[3], [4], [5], [6], [7], [8], [9], [10]. The opportunity to also 
design the excitation input signal to be used in the experiment 
opens up for the possibility to connect the system identifi- 
cation experimental conditions to the required the control 
performance. One way is to formulate this as a convex op- 
timization problem,[ 11], [12], [13], [14], [15], [16]. We will 
study one aspect of the so-called applications oriented input 
design introduced in [17], specifically for model predictive 
control (MPC). The objective is to guarantee, with a given 
probability, that the estimated model belongs to the set of 
models that satisfies the control specifications. This objective 
can be stated mathematically as a set constraint where the set 
of all identified models corresponding to a particular level 
of confidence must lie inside the set of all models fulfilling 
the control specifications [17]. To ensure that the obtained 
optimization problem is convex, we generally must make 
a convex approximation of the set constraint. Two known 
approaches of doing this are the scenario approach, [18] and 
[19], and the ellipsoidal approach, [20]. The main drawback 
of these methods are the computational efforts necessary 
to obtain a descent approximation. Both methods require 
several simulations to be made of the closed loop system 
with MPC. In this paper, we introduce a new method of 
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approximating the set constraint with a convex one. The 
method is based on a perturbation analysis which only 
requires one simulation of the closed loop system. Thus, the 
method is expected to be much faster than both the scenario 
and the ellipsoidal approach. The outline of the paper is 
as follows. In Section [U] we go-through the mathematical 
background necessary. We describe the scenario approach 
and the ellipsoidal approach in Section [III] followed by a de- 



tailed description of the proposed new method in Section IV 
In Section fVl we illustrate the method in two numerical 



examples and in Section VI some conclusions are stated. 



II. Preliminaries 

A. System and model 

We consider a linear, time-invariant, asymptotically stable 
system in discrete time. The system is 



x{k+\) =Ax{k) +Bu{k), 
y(k)=Cx{k). 



(1) 



Here, x € is the state vector, u € W" is the input vector 
and y{k) £ M."y is the output vector The matrices A, B and 
C are the state space matrices of the system. In system 
identification, we want to find a model of the system ([T}. 
We assume that the model is parametrized with an unknown 
parameter vector e M", that is. 



x{k+i,9) =A{e)x{k,e)+B{e)u{k,d), 
y{k,e)^c{e)x{k,e). 



(2) 



In addition, we assume that the model (|2]l matches system ([T]) 
exactly when Q ^ Q,,- We call 0„ the true parameter vector 
The objective of system identification then is to estimate 
the values of 9 that best describes the system, according 
to some measure. The estimated parameter vector, given 
measurements in the identification experiment, is denoted Q^. 

B. Model predictive control 

Model predictive control (MPC), also referred as receding 
horizon control, is an advanced optimization based control 
technique. At each control interval, MPC computes a se- 
quence of optimal inputs by solving an on-line optimization 
problem, where a model is used to predict the behavior of 
the plant. However, only the first input value is applied to 
the plant. A common optimization problem that is solved at 



every time instant t, in MPC is 



min J=J^\\y{k+UO)-r{k+l)\\l+Y,\\Au{k,0)\\l 



{«(<:,e)}™^j k=0 k=\ 

s. t. x{k+\,d) =A{e)x(k,e)+B{e)u{k,e),k^l,...,Ny, 
y{k+l,9)=C{9)x{k+l,e), k = 0,...,Ny, 

x{i,e)=x*it,e), 

Am(1,0) =M(1,0)-M*(f-1,0), 
Umin ^ u[k,9) < M^aj;, A: — 1 , A^„, 
Jmm < 1,0) < ymax,k = 0, ...,A^v 

(3) 

Here r{k) is the reference trajectory, e M" is the vector 
of system parameters, Q and R are weight matrices, A^„ 
and A', are control and prediction horizons, respectively, and 
Au{k,e) = u{k,e)~u{k-l,9). Note that AM(/t,0) = for 
k> Nu- M* (f — 1 , 0) is the optimal input value applied to the 
system at time instant f — 1, and x* (f , 0) is the system state at 
time f, which can be obtained by direct measurement or an 
observer. Different MPC formulations are discussed in more 
detail in [21]. 

C. Prediction error method 

We use the prediction error method (PEM) to estimate the 
unknown parameters of a considered system. The unknown 
parameters are denoted e M", the true parameters repre- 
senting the system are denoted 0^ G M" and the estimated 
parameters given measurements are denoted 0^ G ■ 
A key asymptotic (A^ — ^ oo) property of PEM, is that the 
estimated parameters lie in an identification set with a certain 
probability. The identification set is defined as 



where 7 is a user-defined positive constant. Every parameter 
vector for which the performance degradation is less than 
1 // can be considered as an acceptable parameter from an 
application's point of view. Therefore, we define the set of 
all acceptable parameters, the application set as 



Ssi{a) = <^ : [0 - 0„]'/f (0„)[0 - 0„] < 



(4) 



where the term Xa(n) is the a-percentile of the x^' 
distribution with n degrees of freedom and If is the Fisher 
information matrix. We thus have that 9n G (^siicc) with 
probability a when N ^ °°. For more details, we refer the 
reader to [22]. 

D. Applications oriented input design 

Model-based controllers, such as MPC, use a model in 
order to control a system. Therefore, the control performance 
is affected by any plant-model mismatch. We use the concept 
of an application cost function to relate the plant-model 
mismatch to the performance degradation. We use a scalar 
function of as the application cost and denote it Vapp{0)- 
The cost function is chosen such that its minimum value 
occurs at = 0o. In particular, we assume that Vapp{Q(i) ~ 0. 
Note that if Vapp(0) is twice differentiable in a neighborhood 
of 00, this implies that 

Vappm - , y;,p(0„) - and y„%(0„) > 0, 

see [20]. For a given plant, there is a limit on the maximum 
value of acceptable performance degradation, that is, 

yapp{0) < 7,, (5) 



0(7) = 0|K,pp(0) < 



(6) 



The apphcation set (j6]l has been extensively used in applica- 
tions oriented input design for system identification (see [20], 
[23] and [17]). The main objective of applications oriented 
input design is to provide a tool for designing the input 
signal to be used in the identification experiment such that the 
estimated model guarantees acceptable control performance 
when used in the control design, that is, we want 9^ G 0(7) 
with high probability. This requirement can be formulated 
mathematically as the set constraint 



<^s/(«) c 0(7). 



(7) 



Therefore, the input design problem can be formulated as 
an optimization problem, where (|7]) plays the role of a 
constraint. However, one crucial issue is that while (05/ is 
an ellipsoidal set, the application set can be of any shape. 
Thus, the set constraint (|7]i may not be convex. Two known 
approaches to make a convex approximation of the constraint 
are discussed in the next section. Alternatives to constraint 
(|7]l can be found in [24]. 

III. Application set approximation 

Two methods of approximating the set constraint with a 
convex one are the scenario approach, see [18], [19], and the 
ellipsoidal approach, see [20]. 

In the scenario approach, the application set is described 
by a number, N^, of samples (or scenarios) which are 
randomly chosen from the set. The constraint ^ is then 
replaced by a set of inequalities. 



[0-0o]^/f(0o)[0-0o]> 



N 



Vapp{ek),k=l,...,Nk- (8) 



However, in order to have a good approximation of the appli- 
cation set, the number of samples must be large enough (see 
e.g. [25] for the minimum required number of scenarios). 
This is not easy to satisfy, especially in high dimensional 
and complex plants, since for certain controllers, such as 
MPC, it is not possible to find analytic expressions for Vapp- 
Therefore, a large number of highly time-consuming and 
costly simulations are necessary. 

The ellipsoidal approach is based on a second order Taylor 
expansion of Vapp{9) around 0o, that is. 



Vapp{9) 



^V«pp(0<,)+V„'pp(0)[0-0„] 

-Q.5[e-eofv::pp{e)[o-o,] 

= O + O + O.5[0 - 0,/y" (0)[0 - 0„]. 



(9) 



The apphcation set can thus be approximated by the ellip- 
soidal set 



'^appiY) = 0|[0 - 0./V„„„(0)[0 - 0„] < 



(10) 



The quality of the approximation not only depends on the ap- 
plication cost but also on the value of 7. For sufficiently large 
values of 7, S'app gives an acceptable approximation while 
for smaller values, higher order terms of Taylor expansion 
may need to be considered [26]. However, calculation of the 
Hessian matrix is a challenging task. In many problems it 
is not possible to analytically determine the Hessian of the 
application function due to nonlinearities in the controllers 
that are being used. Therefore, numerical approximations 
are used. Using numerical methods, such as finite difference 
approximation, is not possible in many cases because of the 
large number of variables involved. 

IV. Application set approximation for MPC 

MPC has drawn much attention in control fields, thanks 
to its ability to cope with system limitations. Using MPC, 
we can deal with both input and output constraints explicitly 
during the controller design and implementation. However, 
the resulting explicit solutions for MPC are difficult to deal 
with due to these constraints, which makes it unavoidable 
to use numerical calculations for the approximation of the 
application cost [23]. In this section we present a new 
approach based on analytical methods for the application 
cost approximation for MPC. The proposed approach leads 
to faster estimations of the application sets. 

A. Application Cost Function for MPC 

The application cost function measures the amount of 
performance degradation that stems from plant-model mis- 
match. One reasonable choice of this function for MPC is the 
difference between the measured output when the controller 
is working based on the true parameters, Qq, and when it is 
using perturbed parameters 0, that is. 



M 



Vappie) = - ^ \\y{t, 0„, eo)~yit,e,e„) 



(11) 



t=i 



where M is the number of measurements used, t is time, the 
second argument of y describes the parameters which are 
used by MPC and the third one represents the true system 
parameters [23]. This is shown in Fig.[T] However, in reality 
the true system parameters are not known. Moreover, it is 
not possible to run the process based on perturbed parameters 
and measure the real output, since the plant is then controlled 
using an arbitrary model and it may damage it. Therefore, 
the following approximation of the application cost is used 
[23] 

VappW^^L\\yi'^^^^)-yi''^M\ (12) 

■"^ r=l 

where 9 is the best available estimation of Oo in the linear 
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Fig. 1. The output signal used in the application cost function jl I) . 

approximation of the true system. Thus, the evaluation is 



done in simulations instead of the real plant. In addition, the 
value of 7 in (jSjl is computed based on the idea developed 
in [23] and we explain it later in the result section. 

B. Application Function Approximation 

In order to obtain a convex approximation of the appli- 
cation set, we start by estimating y{t,9,0) in (12i. Using a 
Taylor expansion of y{t, 9, 9), we can write 



yi9)^y{9)- 



1 « " 523^ 



(13) 



where 9j are the elements of 9. Here, the first and third 
arguments of y{t, 9, 9), are omitted for the sake of simplicity. 
In order to find the derivatives in ( [T3| ), we need to find 
the derivatives of the input signal generated by MPC with 
respect to 9. However, this is a challenging problem since 
the solution of MPC is not simple enough when there 
are inequality constraints on input and output signals. The 
proposed solution here is to notice that, when is a small 
perturbation of 9, the active constraints are the same as when 
the MPC is based on 9. Thus, the main idea is to let MPC 
run based on 9 at each time instance t, and determine the 
optimal value of the input signal u{t, 9). We assume that the 
active constraints remain the same for small perturbations 
of 9. Therefore, at time step f, we are able to find an 
explicit solution of the optimization problem in MPC for 
9 = 9 + 59 by considering active constraints as equality 
constraints. We can analyze the effects of perturbing the 
parameters when 50 is small enough. In the rest of this 
section, we briefly describe the explicit solution of MPC 
when we are considering only active constraints, then we 
provide insights into the perturbation analysis for the MPC 
solution. Finally, we show how these concepts can be used 
to find the derivatives in ( [T3] l and compute the application 
cost function. 

1) Explicit Solution of MPC: Consider the MPC problem 
Q at time instance t. For simplicity it is assumed that 
Nu = A^v- Now we seek to rewrite the MPC formulation 
as a quadratic program where we are considering only 
active constraints obtained by solving MPC for 9, which 
are equality constraints. Introducing 

l,9f,...,x{h0f,u{N,,9),...,u{l,9ff, 



X{9) 



-I 





/ / 

/ 



/iV„ + l®C(0) 









In,+i(E>Q 
In„<»R_ 

^=[r(A^„ + l)^,...,r(l)^O,...,O,M*(^-l,0)^]^ 

where by !„ <Xi M, we mean the Kronecker products of Im 
and M [27], we can rewrite the cost function 7 in ([3]) in the 
following form: 

J={Y{9)X{9)-Jf)^{Y{9)X{9)-Jff. (14) 



Moreover, the system dynamics and the first equality con- 
straint in (|3}, give that '^(0)X(0) = ^(0), with 



7-A(0)...O -B(0). 



-m 





(15) 



.../-A(0) 
...0 / 




i(f,0) 

Now consider the inequahty constraints in ([3]l. They can be 
rewritten as 



/iv„+i'^C(0) 

-/^„+i^C(0) 

/ 

-/ 



X(0)< 



In„ €5 Umax 



(16) 



Let E be a diagonal matrix, where each diagonal element 



corresponds to one of the inequality constraints in ( 16 1. A 
diagonal element is zero if its corresponding constraint is 
inactive and it is one for active constraints. Multiplying ( 16 1 
by E and introducing 

" /iv„+i<»C(0) 0^ 
-/w„+i®C(0) 
/ 
-/ 



~^Nu + l ®ymin 



we get Za= p, which represents those inequality constraints 
that are active at time instance t. Then we can rewrite the 
entire set of constraints as J2/(9)X{9) — ^M{9), where 



(0) = 



no) 



'^(0) 
P 



(17) 



Finally, the following optimization problem is obtained: 

min (T(0)X(0)- Jf)^(T(0)X(0)-jr)^, 

s.t. ^(0)X(0) = ^(0). 

Problem ([TTJi is a quadratic optimization problem with equal- 
ity constraints. The KKT conditions [28] for this problem are 

2T(0)^^(T(0)X(0)-^) + ^^(0)A =0, 

£/(0)X(0) =^(0), 

where A are the Lagrange multipliers. This can be written 



as 



2T(0)^^T(0) £/{9y' 
£/{9) 



-\X{9j 




[ ^ 





2Y{9y^Jf{9) 

Sg{9) 



or equivalently 



'P(0) 



X{9] 
A 



= A(0). 



,(18) 



(19) 



Since the block matrices in ^(0) are not invertible, ( 19 1 can 
be solved using the pseudoinverse and Schur complement of 
the resulting block matrix 



X{9) 
A 



(^(0)'^(0))-i^(0)'A(0) 



(20) 



Solving (20 1, we can easily obtain an explicit solution, ^(0), 
for ([TTJl. 

2) Perturbation Analysis: The analysis in this section are 
based on the perturbation analysis techniques in [29] and 
[30]. Having the MPC solution at time step f as a function 
of 0, our aim is to compute the derivatives of X{9) with 
respect to 0, based on which the derivatives in ( [T3| ) will be 
calculated. This can be obtained by linearizing X{9) around 
0, invoking the Taylor expansion 



X(0)=X(0) + f ^^50, 



<90, 

n 



1=1 

1 ^ ^ d^X{9] 
2 



(21) 



d9id9j 



59i59j + hot , 



where = + 50. Moreover, the Taylor expansion of X(0) 
can be computed writing the Taylor expansions of £/{9), 
!3§{9), T(0), and J^(0), which in turn are easily derived 
by having the derivatives of A(0), B{9), C(0), i(f,0), and 
u*{t— 1,0). The derivatives of i(r,0), and u*{t— 1,0) are 
available from the Taylor expansion of X(0) in the previous 
time instances. Now, recall the definition of ^(f , 0, 0) and the 
linear model used for description of the plant 



(22) 



jc(f + 1,0) =A{9)x{t,9)+B{9)u{t,9), 
y{t,9,9)=C{9)xit,9), 

where u{t,9) is the optimal input designed by MPC. We aim 
to find the coefficients in ( pjj l. They can be calculated easily 
in a recursive manner by differentiating (22 1 with respect to 



0, using the derivatives of u{t,9), which are available from 

3) Application Cost Function: Recall the application 
function ( [T2] i, we can calculate the Hessian matrix in terms 
of the obtained derivatives of y as follows 



M 



Z."^ ;in ^ ^ d9 ^ 



=1 



2 ^ 
M 



d9 

,d^y{t,9) 



(23) 



j:{'^f^V{yif^0,9)-y{t,9,9)}. 



Note that the second term is zero since y„,,p(0) = 0. Sub- 
stituting ( [23] ) into ( [TO) l, we a convex approximation of the 
application set. 

The method provides a fast tool for convex approximation 
of application cost function. Many calculations in differ- 
ent time instants are the same and can be pre-computed. 
Moreover, the active constraints may not change often, thus, 
at each time instance a large number of the calculations 
can be skipped by re-using the results from previous time 
instances. Therefore, the proposed approach is much faster 
than both the scenario-based approach and the ellipsoidal 
approximation method. 



V. Numerical examples 
In this section we evaluate the proposed method in Sec- 



tion 



IV with two numerical examples. 



A. Example 1 

Consider the following system: 

x(f + 1) = d2x{t) + u{t), 
y{t) = eix{t). 



(24) 



The true system is given by the parameter values Qq = 
[0.6 0.9]^. The objective is to find the application set 0, 
when MPC is used for reference tracking. We use the MPC 
formulation in (|3]l, with the following settings: A^„ =Ny = 5, 

Q — 10, 7? = 1, Wm«x ~ ^min ^ 1? ymax — ymin ~ 2. 

We set the length of the experiment to = 100 samples 
and the accuracy to 7= 1000. Note that we use the applica- 
tion cost function defined in (111. Now using the proposed 



approach, we obtain the application ellipsoid shown in Fig. 
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Fig. 2. Approximated Eapp ('-') and 400 randomly generated .samples of 
9. '*' represents samples which are satisfying Vnpp(6) < ^ while 'o' are 
those located outside the application set. 8 1 % of the samples that fulfill the 
condition V„pp{Q) < are located inside the approximated ellipsoid given 
by the proposed method. 



In order to check the accuracy of the proposed method, 
we perform 400 simulations with different values of 6 which 
are generated randomly with a uniform distribution. The 
results show that from 400 generated points, 93 points are 
satisfying the condition Vapp{0) < ^. Among all accepted 
values of 9, 81% are completely inside or on the border of 
the approximated ellipsoid, which means that the estimated 
ellipsoid covers at least 81% of the acceptable. 

Furthermore, the Hessian matrix is computed employing 
numerical methods, provided by DERIVESTsuite. The appli- 
cation set is then approximated using the ellipsoidal approach 
(10 1. As expected, the result is the same as when the pro- 



posed method is used. However, in the proposed method, we 
need only one complete simulation of the closed loop system 
with MPC, while in the numerical approximation of the 
Hessian, which is based on finite difference approximation, 
0{6*n^) number of simulation is required depending on the 
selected accuracy. Therefore, the new approach is expected 



to be faster. While it takes 94 seconds for the numerical 
method to calculate the Hessian matrix in this example, 
the new method needs only 12 second to give the same 
approximation, which means that 87% of time is saved. 

B. Example 2 

In this example we illustrate the algorithm on a more com- 
plex and experimental example. We consider a distillation 
column. The nonlinear system representation is taken form a 
benchmark process proposed by the Autoprofit project [31] 
is used. For a general description of distillation columns, we 
refer the reader to [32]. 

The plant is linearized around the steady state operating 
conditions and then, using model order reduction methods, 
the second order model 



x(f + 1) 



01 

03 



02 
04 



xit)- 



05 
07 



0fi 



u{t), 



-0.8954 0.1421 
-0.2118 -0.1360 



(25) 



x(f)+e(f), 



is obtained, where, e(f) is a white measurement noise with 
variance E{e{t)^ e{t)} = 0.001. We assume that 1% perfor- 
mance degradation from the case when MPC is using the 
true parameters is allowed, that is. 



100 



where y(0o)-^I?lilb(f,0o,0o)-Kf)ll^ see [23]. 

Since MPC is used for tracking, the model is augmented 
with a constant output disturbance on each output to get 
integral action. This is presented in further detail in [21]. 
The proposed method has been employed to calculate the 
approximate application cost in (|9]l. In order to evaluate 
the capability of the method, we run the process for 100 
different values of 9, taken from a uniform distribution. Fig. 
[3] shows the real and approximated values of the application 
cost function for each scenario. In order to have a better 
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Fig. 3. Approximated ('o -') and real ('. - -') values of V<,pp(6) for 100 
different samples of 9 taken form a uniform distribution. 

insight, the samples which are located inside the application 
set are illustrated in Fig. |4] It can be easily seen that 
the proposed method has a good performance inside the 



application set. Among 85 scenarios which result in an 
acceptable application cost, 83 scenarios are approximated 
as acceptable ones using the proposed method. The method 
classifies 6 points outside the region as acceptable ones. 
Therefore, the obtained accuracy of the proposed method 
is 92%. 
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Fig. 4. Approximated ('o -') and real ('. - -') values of Vapp{6) inside 
the application set. 92% of the samples inside the region are classified as 
acceptable ones by the proposed method. 



VI. Conclusions 

In this paper we have introduced a general technique for 
the approximation of the application set, a structure required 
for the implementation of optimal input design schemes. In 
particular, we have focused on MPC, a control technique for 
which it is not possible to obtain the application set explicitly. 
Some simulation examples have been presented, which show 
the advantages of the new method with respect to previous 
techniques,in terms of speed. 

The method is general enough to be applied to other 
controller strategies and application areas where it is not 
possible to derive the application set expUcitly. Specifically, 
the method can be extended to MPC for nonlinear plants, 
with more complicated noise structures, and the derivation of 
expressions for higher order derivatives of the cost function 
could be used, in principle, to obtain better approximations of 
the application set using techniques such as the one presented 
in [26]. 
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